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Abstract. In this research note, we provide a variational basis for the optimal artificial diffu- 
sion method, which has been a cornerstone in developing many stabilized methods. The optimal 
artificial diffusion method produces exact nodal solutions when applied to one-dimensional prob- 
lems with constant coefficients and forcing function. We first present a variational principle for 
a multi-dimensional advective-diffusive system, and then derive a new stable weak formulation. 
When applied to one-dimensional problems with constant coefficients and forcing function, this 
resulting weak formulation will be equivalent to the optimal artificial diffusion method. We present 
representative numerical results to corroborate our theoretical findings. 

1. INTRODUCTION 

Many transport-related processes are modeled as advective-diffusive system. For example, trans- 
port of contaminants in subsurface flows is modeled as an advection-diffusion equation coupled with 
Darcy flow. Except for very simple and limited problems, one cannot find analytical solutions and 
hence must resort to numerical solutions. However, it is well-known that great care should be taken 
in developing numerical formulations in order to avoid spurious oscillations due to the advective 
term. 

Many numerical formulations have been proposed, which fall under the realm of stabilized meth- 
ods. One such method is the optimal artificial diffusion method, which has been the basis for 
developing many stabilized methods and has also served as a benchmark for comparison [6l [2]. 
The optimal artificial diffusion method is derived by imposing the condition that it should produce 
exact nodal solutions for one-dimensional problems with constant coefficients and forcing function. 
This derivation, as one can see, does not have a variational basis. 

Herein we outline a variational structure behind the optimal artificial diffusion method. We 
start with a variational statement for the advection-diffusion equation (which is not a self-adjoint 
operator), and derive a stable weak formulation. The resulting weak formulation when applied 
to one-dimensional problems with constant coefficients and forcing function produces the same 
difference equation as the optimal artificial diffusion method, which produces exact nodal solutions 
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when applied to such problems. This shows that the optimal artificial diffusion method has a firm 
variational basis. This paper also highlights other possible routes in developing stable formulations 
for non-self- adjoint operators. 

An outline of this short paper is as follows. In Section [21 we present governing equations for an 
advective-diffusive system, and also outline the standard weighted residual method. In Section[3l we 
describe the optimal artificial diffusion method. In Section [H we briefly discuss Vainberg's theorem, 
which provides a connection between the weighted residual statement and its corresponding scalar 
functional (if it exists) . In Section \5\ we present a variational principle for the advection-diffusion 
equation, and then derive a stable weak formulation. We then show that the resulting weak 
formulation produces the same difference equation as the optimal artificial diffusion method for one- 
dimensional problems with constant coefficients and forcing function. Finally, we draw conclusions 
in Section [6l 

Remark 1.1. It should be noted that there is a huge literature on developing stabilized (finite 
element) formulations for an advective-diffusive system. For example, see ^''^'^ references 

therein. A thorough discussion of these works is beyond the scope of this paper. In addition, these 
discussions are not relevant to the subject matter of this paper as none of them discuss a variational 
principle (that is, constructing a scalar functional) for an advection-diffusion system. 

2. GOVERNING EQUATIONS: ADVECTIVE-DIFFUSIVE SYSTEM 

Let C M"^ be a smooth and bounded domain, where "d" denotes the number of spatial di- 
mensions, and dft denotes its smooth boundary. As usual, the boundary is divided into F^ (the 
part of the boundary on which Dirichlet boundary condition is prescribed) and F^ (the part of 
the boundary on which Neumann boundary condition is prescribed) such that F^ U F^ = dil and 
-pD PI -pN _ 0^ -pgi. ^ J7 denote the position vector, and the gradient and divergence operators 
are denoted as "grad" and "div", respectively. Let u : 17 — > M denote the concentration, v(x) the 
velocity vector field, and k{x) the symmetric and positive-definite diffusivity tensor. For further 
discussion consider the following steady advective-diffusive system 

(2.1) v{x) ■ grad[u(a;)] — div [k{x) grad[u(a;)]] = f{x) in Q 

(2.2) u{x) = u'P{x) onF° 

(2.3) n{x) ■ k{x)giad[u{x)] = t^{x) on F^ 

where u'^{x) is the prescribed Dirichlet boundary condition, t^{x) is the prescribed Neumann 
boundary condition, f{x) is the prescribed volumetric source, and n{x) denotes the unit outward 
normal to the boundary. Note that the advective-diffusive operator is not self-adjoint. It is well- 
known that, even under smooth functions for k{x),v{x) and f{x), the solution u{x) to equations 
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(j2.ip ~ (j2.3p may exhibit steep gradients close to the boundary especially for advection-dominated 
problems [B] (and also see Figure [TJ which will be described later in Section [3]). 

2.1. Notation and preliminaries. In the next section, we present a weak formulation for an 
advective-diffusive system. To this end, let us define the following function spaces: 

(2.4a) U := {u{x) G H^{n) \ u{x) = uP{x) on r°} 

(2.4b) W := {w{x) e H^{n) \ w{x) = on r°} 

where H^(Q) is a standard Sobolev space defined on Q. Note that the inner product for the above 
vector spaces is the standard inner product. That is. 



(2.5) (a; 6) := / a{x) ■ h{x) m 

Jn 

Similarly, one can define the weighted inner product 

(2.6) {a;b)^:= / fi{x)a{x) ■ b{x) dQ 



where fi : Q ^ M"^ (where denotes the set of positive real numbers) is the (scalar) weight 
function or measure density. Note that a weight function (which is used to define the weighted 
inner product) should not be confused with weighting functions (which are sometimes referred to 
as test functions) . Using this weighted inner product one can define weighted Sobolev spaces as 

(2.7a) Llin) := [u{x) \ (n; n)^, < +00} 

(2.7b) Hlin) := [u{x) G Ll{n) \ (grad^; grad^)^ < +00} 

Corresponding to the function spaces given in (j2.4p . we can define the following weighted function 
spaces 

(2.8a) := {u{x) G Hj^{n) \ u{x) = uP{x) on r°} 

(2.8b) W^, := {w{x) G Hj^{n) \ w{x) = on T^} 

2.2. Standard weighted residual method and the Galerkin formulation. Let w{x) denote 
the weighting function corresponding to u{x). A weak formulation based on the standard weighted 
residual method for the advective-diffusive system given by equations (|2.ip - (j2.3p can be written as 



(2.9) Find u{x) G U such that T{w; u) = yw{x) G W 

where the bilinear functional J-" is defined as 

T{w;u) := / w{x)v{x) ■ gTad[u{x)] + / grad [w(a;)] • fc(a;)grad [n(a;)] dfl 
Jn Jn 

(2.10) - / w{x)f{x) dn- [ w{x)tT'{x) dr 

Jn JrN 



A finite element formulation corresponding to the standard weighted residual formulation can be 
written as 

(2.11) Find u'^(a;) G such that T{w^]u^) = 'iw^{x) G W'' 

where lA^ and W'* are finite dimensional subspaces of lA and W, respectively. One obtains the 
Galerkin formulation by using the same function space for both hi^ and (except on F^). 

The Galerkin formulation produces node-to-node spurious oscillations for advection-dominated 
problems [6]. The Galerkin method loses its best approximation property when the non-symmetric 
convective operator in the transport equation: for example, see the performance of the Galerkin 
formulation in Figured] In other words, the Galerkin method is not optimal for solving advection- 
dominated problems. In principle, it is possible to choose a small enough grid such that the 
element Peclet number is less than one, and avoid spurious oscillations under the Galerkin method. 
However, it may not always be practical to choose such a fine grid, and therefore, one needs to 
employ a stabilized formulation to avoid unphysical oscillations and get meaningful results on coarse 
grids. To understand this anomalous behavior several theoretical and numerical studies have been 
performed, see [H O [2] and references therein. One of those studies is a simple method that gives 
nodally exact solutions for one-dimensional problems with constant coefficients, which is commonly 
referred to as the optimal diffusion method. 



3. A NODALLY EXACT FORMULATION AND OPTIMAL ARTIFICIAL DIFFUSION 

Consider the following one-dimensional advection-diffusion equation with homogeneous Dirichlet 
boundary conditions, and constant coefficients (that is, velocity, diffusivity and forcing function are 
constants) : 

du , cfiu 

(3.1) -T.-^lUS-l ^"("•^' 

(3.2) u{x = 0) = and u{x = 1) = 

The analytical solution for the above problem is given by 

/o Q^ I \ f ( ^~ exp(ua;/fc) 

(3.3) u{x) = — \x 



V \ 1 — ey.-piv/k) 

which is plotted in Figure [T] for various value of v/k. As one can see, for large v/k we have 
steep gradients near the outflow boundary. The above problem has been used as a benchmark for 
developing many stabilized formulations [B|. 

For a numerical solution, let us divide the unit interval into N equal-sized elements (and hence 
+ 1 nodes), and define h := 1/A^. Let the nodes be numbered as j = 0, • • • , A^. Then, the position 



vector of node j is xj = jh. The difference equation at an intermediate node {j = 1, • • • , — 1) 
arising from the Galerkin formulation for the above problem can be written as 

(3.4) ^ l"^"-^ + p-f^ + j = / 

where = vh/{2k) is the element Peclet number, and uj is an approximate numerical solution at 
node j. That is, 

(3.5) Uj ^ u{xj) Vj = 0, • • • , 



The difference equation ()3.4p can be rearranged as 

(3.bj V - fc p - / 

The above equation basically reveals that the Galerkin formulation approximates the first- and 
second-derivatives using a central difference approximation at intermediate nodes. 

It is well-known that the Galerkin formulation is under-diffusive when applied to the advection- 
diffusion equation, which is considered to be the reason why the formulation gives node-to-node 
spurious oscillations. In Figure [H we compare the numerical solution from the Galerkin formulation 
with the analytical solution for the aforementioned one-dimensional problem. As one can see from 
the figure, the Galerkin formulation produces spurious node-to-node oscillations for high Peclet 
numbers. 

To also understand this anomalous behavior of the Galerkin formulation, we now outline a nodally 
exact formulation for the above problem, which is commonly referred to as the optimal artificial 
diffusion formulation in the literature. To this end, we start with a difference equation at node j 
of the form 

(3.7) + PoUj + (3iUj+i = f 

and the unknown coefficients (/3_i, /3o, and /3i) are determined so that the above difference equation 
gives nodally exact solution for the model problem given by equations (|3.ip - (|3.2p on the uniform 
computational mesh described above. After simplification, the difference equation for the nodally 
exact formulation can be written as: (see References [?l \6\ \TU\), 



(3.8) 



V 

2h 



[1 + cothPg^)Mj_i + {2cothP^)uj + coth Pl^)uj 



+1 



/ 



where "coth" denotes the hyperbolic cotangent function. By rearranging the terms one can write 
the above equation as follows 



(3.9) 



2h 



where the artificial diffusion coefficient k (which depends on the mesh size, and medium and ffow 
properties) is given by 

(3.10) k:=^(^cothP^-^^ 

The difference equation (j3.9p reveals that the nodally exact formulation (or the optimal artificial dif- 
fusion method) also employs a central difference approximation for the first- and second-derivatives 
but solves a modified advective-diffusive system with an additional artificial diffusion given by the 
coefficient k. 

It should be noted that the above derivation for the optimal artificial diffusion method is not 
based on a variational principle. In this paper we present a variational structure behind the 
optimal artificial diffusion method, which, to the best of our knowledge, has not been reported in 
the literature. We start by presenting a variational principle for an advective-diffusive system and 
then derive a weak formulation. A finite element approximation of this weak form gives the optimal 
artificial diffusion method. 

We now present a variational principle for an advective-diffusive system. We then outline the 
underlying variational structure behind the optimal artificial diffusion method. 



4. VAINBERG'S THEOREM AND EXISTENCE OF A SCALAR FUNCTIONAL 



Vainberg's theorem [K] provides a connection between scalar functionals (also referred to as 
"energy" functionals) and weighted residual statements. The theorem provides a criterion to es- 
tablish when a scalar functional exists. The theorem also provides a formula to compute the scalar 
functional (if exists) from the weighted residual statement, which we will not invoke in this paper. 

Vainberg's theorem can be stated as follows. Let Q{w;u) be a weighted residual functional, 
which is linear with respect to w (but need not be linear with respect to u). There exists a scalar 
functional 8{u) such that 



Q{w;u) 



Eiu + ew) 

de 



e=0 



if and only if 
(4.1) 



de 



G{w, u 



ew 



e=0 



de 



G{w, u + ew) 



e=0 



For a simple proof of Vainberg's theorem see Reference [9l Chapter 9]. 

Returning back to our discussion on advective-diffusive system, using Vainberg's theorem one 
can conclude that there exists no scalar functional £{u) whose directional derivative along w{x) 
gives the (bilinear) functional in the standard weighted residual method (which is given by equation 



(j2.10p ). That is, there is no scalar functional £{u) such that 



J^iw; u) 



— £{u + ew) 



J e=0 



However, we should note conclude that the advective-diffusive system does not have a variational 
principle. All we concluded above (invoking Vainberg's theorem) is that the standard weighted 
residual statement (j2.10p cannot be obtained as a directional derivative of a scalar functional. It 
may be possible to construct a different weighted residual statement with different weight function 
(or measure) that can be obtained as a directional derivative of a scalar functional. Stated differ- 
ently, just because the operator is non-self-adjoint does not mean that there exists no variational 
statement. 

In a seminal paper, Tonti |14j has correctly highlighted the point that the symmetry of a bilinear 
form in the sense of equation (j4.ip (which according to Vainberg's theorem is a necessary and 
sufficient condition for the existence of a variational formulation) depends on the choice of the inner 
product. (Note that Tonti has used the acronym "bilinear functional" to mean inner product.) One 
of the ways to meet the symmetry requirement is by changing the underlying inner product. Tonti 
|14j has also shown that changing the inner product is equivalent to transforming the given problem 
into a different problem that has the same solution(s). Another related work is by Thangaraj and 
Venkatarangan [TH] who have presented dual variational principles with applications to magneto- 
hydrodynamics. A related and important work is by Magri [11] who has developed a procedure to 
select an appropriate inner product for any given linear operator to meet the symmetry requirement. 
Following the discussion in References [14^ [13] , in the next section we present a variational principle 
for an advection-diffusion system. 

Remark 4.1. For self-adjoint operators (e.g., the Laplacian operator) one can show that the stan- 
dard weighted residual method can be obtained as a directional derivative of a scalar functional. 

5. A VARIATIONAL PRINCIPLE FOR THE ADVECTION-DIFFUSION EQUATION 
Consider the following constrained extremization problem 



"(a:) Jn 

where F : 0, X M X M'^ ^ M is a known scalar functional. It is well-known (for example, see [5} 110]) 
that the Euler-Lagrange equation for the above extremization problem is given by 



where F^ denotes the partial derivative with respect to u{x), and Fp denotes the vector of partial 
derivatives with respect to the components of grad[n(a;)]. 

7 



(5.1) 




(5.2) 



Fu {x, u{x),gra.d[u{x)]) — div [Fp{x, u{x),gra.d[u{x)])] = Va; G 



We now construct a scalar functional for the advection-diffusion equation. To this end, let us 
define the scalar function a{x) as 



(5.3) a{x) := exp 



(fc ^{x)v{x)^ ■ dx 



Note that the integral inside the exponential operator in equation (|5.3|) is an indefinite line integral. 
Since the tensor k{x) is positive-definite, its inverse exists, and hence the scalar function a{x) is 
well-defined. It is important to note that a{x) > Va; € Q (as we have assumed the domain to be 
bounded and smooth). Now define the functional for the advection-diffusion equation as 

(5.4) F(a;, n(a;), grad[n(a;)]) := a{x) ^\\y^k{x)gTad [u{x)] ||^ — u{x)f{x)^ Va; G 

where || • || denotes the standard Euclidean norm, and y^k(x) denotes the square root of the 
symmetric and positive-definite tensor k{x). (Note that the square root of a symmetric positive- 
definite tensor is well-defined, and is also symmetric and positive-definite. See Halmos [8] and also 
Gurtin [71 page 13].) A straightforward calculation shows that the Euler-Lagrange equation of the 
functional (15.41) is 



(5.5) v{x) ■ grad[u(a;)] — div [k{x) grad[ti(a;)]] = f{x) 

which is the (multi-dimensional) advection-diffusion equation. For the advective-diffusive system 
given by equations (|2.ip - (j2.3p (that is, including the boundary conditions), the variational statement 
can be written as 

(5.6) extremize I{u) 

u{x) 

(5.7) subject to u{x) = u'^{x) on r° 
where the functional I{u) is defined as 

(5.8) l{u):= / F{x,u{x),gTad[u{x)]) dn- / a{x)u{x)t^{x) dV 

and the scalar function a{x) is same as before, equation (|5.3p . 

We now use this variational principle to derive a new weak formulation for the advective-diffusive 
system. We later show that for constant coefficients and homogeneous Dirichlet boundary condi- 
tions, this formulation is same as the optimal artificial diffusion method. 

Remark 5.1. It is well-known that, using the least-squares method one can construct a minimiza- 
tion problem for a given partial differential equation of the form Cu — f = as follows: 



(5.9) minimize 

u{x) 



[ \\Cu-ffdn 

Jn 



But it should be noted that the Euler- Lagrange equation of the above minimization problem need 
not be the original differential equation Cu — / = 0. For example, if we consider the differential 
equation to be d^u/dx^ = in Q = (0,1), the corresponding minimization problem based on the 
least-squares method takes the following form: 

fd^u\^ 

(5.10) minimize j ( ^^-^ J dx 

However, the Euler- Lagrange equation of the above minimization problem is d^u/dx^ = (which 
is not same as the differential equation we considered). Also, another difference between the mini- 
mization problem based on the least-squares method and the variational problem is in the regularity 
requirements. For example, in the minimization problem based on the least-squares method (|5.10p . 
u should be twice differentiable (or, more precisely, u G -ff^((0, 1))J. On the other hand, the (stan- 
dard) variational principle for the differential equation (f'u/dx^ = requires that u be differentiable 
once (that is, u G LI^[{0, 1)) ). 

5.1. A new stable formulation. A necessary condition that the extremum of the functional X(u) 
has to satisfy is 



(5.11) 



^ I{u{x) + ew{x)) 



= yw{x) 

e=0 



For convenience let us define the bilinear functional B{w; u) as 



B{w; u) := 



I{u{x) + euj{x)) 

de 



=0 



a{x)grad[w{x)] ■ (fe(a;)grad [n(a;)]) dO 



(5.12) - / a{x)w{x)f{x) dn- [ a{x)w{x)tP{x) dT 
A new weak formulation for the advective-diffusive system can be written as 

(5.13) Find u{x) G Ua such that B{w; u) = yw{x) G >V„ 



where Ua and are weighted function spaces defined in equation (j2.8p with weight function 
fi{x) = a{x). Recall that the scalar function a{x) is defined in equation (jS.lip . A corresponding 
finite element formulation can be written as 

(5.14) Find u'^ix) G Uj^ such that B{w^;u'') = Vw''(a;) G 

where and are finite dimensional subspaces oiUa and Wo, respectively. It is important to 
note the formulation ()5.13p (and hence the formulation ()5.14p ) is valid even for spatial dimensions 
nd= 1,2,3. 



Remark 5.2. Some notable differences between the bilinear functions J- {w,u) and B{w;u) are as 
follows. A non- symmetric term similar to w{x)v{x) •grad[n(a;)] d$7 is not present in the bilinear 
functional B{w;u). The bilinear functional B{w,u) is symmetric in the sense that 



— B{w, u + ew) 



— B{w; u + ew) 



e=0 



which is not the case for the bilinear functional T{w]u). The symmetry ofB{w;u) (which, according 
to Vainberg's theorem, is necessary and sufficient for the existence of a scalar functional) should 
not be surprising as we defined the bilinear functional B{w, u) (j5.12p as a directional derivative of 
a scalar functional. Another notable difference is that the bilinear functional J-{w; u) is based on 
the standard inner product, which has the weight function (or measure density) to be 1. The 
bilinear functional B{w; u) is based on the weighted inner product with the weight function equal 
to a{x) > (which is defined in equation (|5.3p ]. 



5.2. Application to the one-dimensional problem. We now apply the new formulation to the 
one-dimensional problem defined in Section [31 and compare the difference equation produced by 
this new stable formulation with the one produced by the optimal artificial diffusion method given 
by equation (|3.9|) for an intermediate node. For the case at hand, the coefficients v and k are 
constants, and hence a{x) = exp[—vx/k]. The variable u{x) is interpolated as 

(5.15) u{x) = Uj-iN^i{x) + ujNq^x) + Uj+iN^i{x) xj-i < x < xj+i 



where the shape functions are defined as 
(5.16a) N_i{x) = I 

(5.16b) No{x) = I 

(5.16c) iV+i(x) = I 



< x - Xj < +h 
^^Zlll+Il -h<x-Xj <0 



h—{x—Xj) 



{) < X — Xj < +h 



h w \ J. 



-h<x-Xj <0 
< X - Xj < +h 



The weighting function is also similarly interpolated. The difference difference equation produced 
by the new stable formulation at an intermediate node can be written as 



(5.17) 
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The coefficients 7_i, 70 and 74.1 can be written as 

J^ll; a{x)k{x)N'Q{x)N'_-^{x) dx 

(5.18a) 7-1 = 

a{x)NQ{x) dx 

j^l:a{x)k{x) K(x))' dx 
(5-18b) 70 = -j^ ^ ^ 

J-h Oi{x)NQ{x) dx 

f^h a{x)k{x)N^{x)N'_^^{x) dx 



(5.18c) 7+1 



J^l^ a{x)No{x) dx 

where a superscript prime denotes a derivative with respect to x. After simphfication, the coeffi- 
cients can be written as 

(5.19) j_, = -^(l + cotHP^)) ,10 = 1 coth(Pi^), 7+1 = ^{l- coth(Pi^) 

As one can see, the coefficients 7_i, 70 and 7+1 are, respectively, the same as the coefficients 
/3o and /J+i (see equations (|3.7p and (j3.8p ). which are obtained using the optimal artificial 
diffusion method. In Figure [21 we compare the numerical solutions obtained using the new stable 
formulation with the analytical solutions. As predicted by the theory, the new stable method 
produces nodally exact solution for all Peclet numbers. Both the theoretical and numerical studies 
have shown that the stable weak formulation proposed in the previous section is equivalent to the 
optimal artificial diffusion method. Hence, we have provided a variational basis for the artificial 
diffusion method, which has academic importance, and also has been the basis in developing and 
testing many stabilized formulations. 

Although the new stable formulation is valid even in two- and three-dimensions, preliminary 
numerical studies have indicated that this formulation may not be practically feasible for large- 
scale problems due to the following reasons: evaluating the exponential function is expensive, one 
may require more Gauss points per element in order to evaluate integrals in the weak form, and 
also the resulting stiffness matrix will be ill-conditioned without an appropriate preconditioner. 
However, the new stable formulation, as shown in this paper, does have theoretical significance. It 
also highlights a possible route of developing robust stable formulations. This paper also highlights 
the relevance of the theoretical studies by Tonti on variational principles (for example. Reference 
|14j ) to develop new stabilized formulations. 

6. CONCLUSIONS 

It is well-known that the classical Galerkin formulation (which is based on the standard weighted 

residual method) produces spurious node-to-node oscillations for the advection-diffusion equation 

for advection-dominated problems. For self-adjoint operators (e.g., the Laplacian equation), the 
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standard weighted residual formulation can be obtained as a directional derivative of a scalar 
functional. However, for the advection-diffusion equation (using Vainberg's theorem) it can be 
shown that the standard weighted residual method cannot be obtained as a directional derivative 
of a scalar functional. 

In this paper, we presented a variational principle for an advective-diffusive system, and de- 
rived a stable weak formulation. This resulting weak formulation when applied to one-dimensional 
problems gives rise to the same difference equation as the optimal artificial diffusion method, which 
produces exact nodal solutions when applied to one-dimensional problems with constant coefficients 
and forcing function. Hence, we have provided a variational basis for the optimal artificial diffusion 
method, which has been the cornerstone in developing many stabilized methods. We presented 
representative numerical results to corroborate our theoretical findings 
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Figure 1. In this figure, we compare the Galerkin formulation (which is denoted 
using solid squares and dotted lines) with the analytical solution (which is denoted 
using solid continuous lines) for one-dimensional advection-difFusion equation with 
homogeneous Dirichlet boundary conditions for various v/k ratios. We have taken 
the forcing function to be unity (i.e., / = 1), and h = 0.1. As expected, the Galerkin 
formulation produced spurious node-to-node oscillations for high Peclet numbers. 
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Figure 2. In this figure, we compare the numerical solution using the new sta- 
ble formulation (which is denoted using solid squares) with the analytical solu- 
tion (which is denoted using solid continuous lines) for one-dimensional advection- 
diffusion equation with homogeneous Dirichlet boundary conditions for various v/k 
ratios. We have taken the forcing function to be unity (i.e., / = 1), and /i = 0.1. As 
predicted by the theory, the new stable formulation produces nodally exact solutions 
for the chosen one-dimensional problem for all Peclet numbers. The continuous line 
denotes the analytical solution, and solid squares denote the numerical solution from 
the new stable formulation. 
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